# Matrix eQTL by Andrey A. Shabalin
# http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/
# 
# Be sure to use an up to date version of R (Revolution R and Matlab are faster on PC).
#
# Set working directory:
# setwd('');
#
library(methods);
args <- commandArgs(TRUE)
setwd(args[1])
source('/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/Matrix_eQTL_R_Source/Matrix_eQTL_R/Matrix_eQTL_engine.R');
SNP_file_name ="SNP.txt" 
snps_location_file_name="/data2/bsi/RandD/bora/OV_EQTL/snpsloc.txt"
# Gene expression file name
expression_file_name = "GE.txt";
gene_location_file_name = '/data2/bsi/RandD/bora/OV_EQTL/geneloc.txt';

# Covariates file name
# Set to character() for no covariates
#covariates_file_name = 'gCovariates.txt';
covariates_file_name =character();
# Output file name
output_file_name = "eQTL_results_R_cis.txt";

# Only associations significant at this level will be saved
pvOutputThreshold_cis = 1e-4;
pvOutputThreshold_tra = 1e-4;

# Error covariance matrix
# Set to character() for identity.
errorCovariance = character();
# errorCovariance = read.table('gerrorCovariance.txt');

cisDist = 1e6;
useModel = modelLINEAR; # modelANOVA or modelLINEAR
## Load genotype data

snps = SlicedData$new();
snps$fileDelimiter = '\t'; # the TAB character
snps$fileOmitCharacters = 'NA'; # denote missing values;
snps$fileSkipRows = 1; # one row of column labels
snps$fileSkipColumns = 1; # one column of row labels
snps$fileSliceSize = 2000; # read file in pieces of 2,000 rows
snps$LoadFile(SNP_file_name);

## Load gene expression data

gene = SlicedData$new();
gene$fileDelimiter = '\t'; # the TAB character
gene$fileOmitCharacters = 'NA'; # denote missing values;
gene$fileSkipRows = 1; # one row of column labels
gene$fileSkipColumns = 1; # one column of row labels
gene$fileSliceSize = 2000; # read file in pieces of 2,000 rows
gene$LoadFile(expression_file_name);

## Load covariates

cvrt = SlicedData$new();
cvrt$fileDelimiter = '\t'; # the TAB character
cvrt$fileOmitCharacters = 'NA'; # denote missing values;
cvrt$fileSkipRows = 1; # one row of column labels
cvrt$fileSkipColumns = 1; # one column of row labels
cvrt$fileSliceSize = 2000; # read file in one piece
if(length(covariates_file_name)>0) {
	cvrt$LoadFile(covariates_file_name);
}
#cvrt=character()
## Run the analysis
snpspos = read.table(snps_location_file_name, header = TRUE, stringsAsFactors = FALSE);
genepos = read.table(gene_location_file_name, header = TRUE, stringsAsFactors = FALSE);

Matrix_eQTL_engine_cis(snps,gene, cvrt, output_file_name,pvOutputThreshold_cis,pvOutputThreshold_tra,useModel,errorCovariance,verbose = TRUE,snpspos,genepos,cisDist)
data<-read.table(output_file_name,sep="\t",head=T)
#s5<-table(data[,6])
#write.table(s5,"out.txt",sep="\t",col.names=FALSE,row.names=FALSE)
data<-data[with(data, order(FDR)), ]
write.table(data,"FDR.txt",sep="\t",col.names=TRUE,row.names=FALSE)
data<-data[with(data, order(p.value)), ]
write.table(data,"PVAL.txt",sep="\t",col.names=TRUE,row.names=FALSE)
